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Abstract:  Recordingsof  seismic  wavefieldsfrom  various 
sources  were  obtained  using  a  small  array  of  vertical- 
component  geophones  under  winter  conditions  at 
Grayling,  Michigan.  These  data  were  processed  using  a 
frequency-wavenumber  domain  Capon  minimum 
variance  beamformerto  estimate  the  bearing  angle  and 
propagation  velocity  of  the  waves  emitted  from  the 
source.  Thecross  power  matrixwas  adaptively  estimated 
using  a  tapered  block-averaging  procedure.  The  wave 
sources  were  sledgehammer  blows  on  the  ground 
surface,  .45  caliber  blank  pistol  shots,  and  an  M60tank 
moving  at  4.5  m  s-1  along  a  road  near  the  array. 
Reliable  wavenumber  spectra  were  obtained  for  all 
sources.  Processing  results  forthe  hammer  blows  show 
that  the  dominant  seismic  arrival  is  a  Rayleigh  wave 
traveling  at  roughly  220  m  s-1 .  Forthe  pistol  shots,  two 
arrivals  corresponding  to  the  airwave  (338  m  s-i)  and 
the  air-coupled  Rayleigh  waves  (220  m  s-i)  were 


observed.  The  results  tor  these  sources  were  relatively 
insensitive  to  the  processing  parameters  used.  For  the 
moving  vehicle,  the  dominant  signals  observed  were 
Rayleigh  waves  (220  m  s-1).  Accurate  locations  were 
obtained  forthis  moving  source,  although  the  processing 
parameters  had  to  be  carefully  selected,  and  the  choice 
of  frequency  parameters  affected  the  accuracy  of  the 
wavenumber  results.  Maximizing  the  number  of  degrees 
of  freedom  and  the  coherence  of  the  frequency  estimates 
and  minimizing  the  variation  of  the  coherence  across 
adjacent  frequency  bins  provided  the  most  consistently 
reliable  strategy  for  obtaining  accurate  wavenumber 
estimates  for  the  moving  vehicle.  The  sensitivity  of  the 
wavenumber  estimates  to  the  frequency  processing 
parameters  seems  to  be  related  to  the  bias  in  the  phase 
spectra  of  the  signals  and  will  potentially  occur  in  any 
bearing  estimation  method  that  uses  temporal  frequncy 
phase  spectra. 
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Source  Location  and 
Tracking  Capability  of  a 
Small  Seismic  Array 

MARK  L.  MORAN  AND  DONALD  G.  ALBERT 


1.  INTRODUCTION 

Acoustic  and  seismic  sensor  array  processing 
has  many  potential  applications  of  interest  to  the 
Army.  Recent  military  conflicts  have  clearly 
demonstrated  the  efficiency  of  smart  weapons 
systems.  A  critical  element  in  the  effectiveness  of 
many  of  these  systems  is  their  ability  to  discern 
and  track  a  target  in  the  presence  of  background 
noise.  The  Wide  Area  Mine  (WAM)  being  devel¬ 
oped  by  the  U.S.  Army  is  heavily  dependent  upon 
acoustic  and  seismic  sensor  information  for  bear¬ 
ing  determination  of  targets.  Weapon  systems  such 
as  the  WAM  are  designed  to  recognize,  track,  and 
destroy  hostile  military  vehicles  such  as  armor, 
mobile  artillery,  and  heavy  transport  vehicles  by 
deploying  small  arrays  of  passive  acoustic  and/or 
seismic  sensors  (microphones  and  geophones)  to 
monitor  local  seismic  and  acoustic  wavefields.  A 
cartoon  of  the  detection-attack  scenario  for  the 
WAM  device  is  shown  in  Figure  1. 

Other  important  areas  where  acoustic  and  seis¬ 
mic  array  processing  may  be  applied  include  pe¬ 
rimeter  defense  and  intrusion  detection.  In  many 
instances  the  effectiveness  of  line-of-sight  systems 
such  as  radar,  infrared,  and  motion  detection  sys¬ 
tems  are  compromised  or  defeated  by  local  topog¬ 
raphy,  vegetation,  and  obscuring  meteorological 
conditions.  In  real-world  settings,  acoustic  and 
seismic  waves  readily  propagate  along  bending 
ray  paths  and  are  less  affected  by  vegetation  and 
airborne  obscurants.  This  allows  coverage  of  zones 
currently  outside  the  detection  limits  of  systems 
that  rely  on  line-of-sight  propagation  paths.  Acous¬ 
tic  and  seismic  array-based  systems  are  relatively 
inexpensive  and  may  be  produced  and  deployed 
in  large  numbers,  providing  a  substantial  perim¬ 
eter  defense  or  detection  capability. 

For  small  array  processing  systems  to  function 
reliably,  it  is  necessary  to  be  able  to  estimate  accu¬ 
rately  the  wavefield's  temporal  frequency  phase 
spectrum,  even  with  limited  spatial  and  temporal 
wavefield  samples.  When  using  wavenumber  es¬ 
timation  procedures  based  on  beamforming  prin¬ 


ciples,  these  data  constraints  make  it  difficult  to 
obtain  reliable  estimates  of  the  spatial  correlation 
matrix,  which  is  critical  to  the  beamformer's  per¬ 
formance.  The  spectral  bias  and  variance  are  the 
statistical  error  parameters  that  describe  the  reli¬ 
ability  of  the  resulting  estimated  frequency  and 
wavenumber  spectra.  Bias  is  a  measure  of  the 
energy  leakage  from  adjacent  frequency  bins  that 
erroneously  appears  at  an  incorrect  frequency; 
this  phenomenon  is  termed  leakage  and  is  caused 
by  the  sampling  distribution  in  the  time  (or  space) 
domain.  Variance  is  a  measure  of  the  stability  in  the 
estimate  and  is  a  function  of  the  number  and 
degree  of  independence  between  samples  avail¬ 
able  to  the  estimation  process. 

Beamforming  complications  may  result  from 
source  motion,  acoustic-to-seismic  coupling, 
inhomogeneous  propagation  environments,  mul¬ 
tiple  sources,  and  low  signal-to-noise  ratios  (SNRs). 
Given  the  nature  of  the  moving  seismic-acoustic 
sources  being  considered  and  the  limited  quantity 
of  data  available,  it  is  important  to  recognize  the 
nonideality  of  the  observed  signals,  so  that  those 
properties  of  the  signal  that  have  the  highest  de¬ 
gree  of  coherency  across  the  observation  array  can 
be  exploited  fully. 

The  combined  effect  of  these  source-dependent 
characteristics  and  the  problems  arising  from 
sparsely  sampled  short-duration  signal  vectors 
has  not  been  previously  demonstrated  on  the  beam 
response  function.  In  this  report  we  show  that 
beamforming  techniques  can  be  reliably  applied 
to  these  types  of  signal  vectors  as  long  as  careful 
attention  is  paid  to  the  spectral  estimation  param¬ 
eters  used  in  the  Fourier  transform  process.  We 
present  a  few  rough  guidelines  that  produce  con¬ 
sistently  reliable  wavenumber  estimates  from  Ca¬ 
pon  maximum-likelihood  beamformers. 

Historically,  frequency  wavenumber  domain 
(F-K)  beamforming  analysis  has  been  applied  to  a 
variety  of  problems.  In  seismology,  wavenumber 
beamformers  using  geophone  arrays  are  routinely 
used  to  determine  source  parameters  for  regional 
and  teleseismic  events  as  well  as  to  identify  nuclear 


explosions.  The  application  of  the  Capon  maxi¬ 
mum-likelihood  beamf  ormer  to  teleseismic  events 
using  the  Large  Aperture  Seismic  Array  in  Mon¬ 
tana  was  first  reported  by  Capon  et  al.  (1967).  In 
1969  Capon  derived  the  beamformer  and  noise 
suppression  characteristics.  Beamforming  is  also 
applied  to  sonar  problems  where  passive  arrays  of 
hydrophones  locate  and  track  submarines  and 
surface  vehicles  (Hahn  1975).  Aircraft  tracking 
with  both  acoustic  and  radar  technologies  and 
stellar  mapping  of  distant  radio-wave  emissions 
are  other  prominent  areas  of  application  for 
beamforming  technology. 

The  use  of  seismic  and  acoustic  arrays  to  track 
and  locate  targets  of  military  interest  dates  from 
the  Vietnam  War.  Beginning  in  the  late  1960s, 
military  efforts  focused,  with  limited  success,  on 
locating  artillery  positions.  In  the  later  1970s, 
targets  such  as  surface  vehicles  were  investi¬ 
gated  (Cress  1976).  Recently,  bearing  estimation  of 
slow-moving  military  aircraft  using  acoustic  and 
seismic  array  processing  technologies  has  been 
actively  investigated.  Of  particular  note  are  the 
efforts  of  Nawab  et  al.  (1985)  and  Lacoss  et  al. 
(1991),  who  have  successfully  tracked  both  air¬ 
borne  and  ground  targets  at  distances  of  a  few 
kilometers. 

Several  researchers  have  investigated  the  im¬ 
pact  of  nonideal  source  attributes  on  bearing  esti¬ 
mates.  De  Graaf  and  Johnson  (1985)  discuss  the 
bias  of  Bartlett  and  Capon  maximum-likelihood 


beamformers  as  a  function  of  SNR  under  the  as¬ 
sumption  that  the  spatial  correlation  matrix  is 
known  exactly  and  the  sources  are  ideal  stationary 
plane  waves.  For  two  or  more  sources  they  show 
that  the  bias  in  the  wavenumber  estimate  is  a 
function  of  the  source  separation  and  the  degree  of 
coherence  between  the  sources,  the  SNR,  the  rela¬ 
tive  source  power,  and  the  number  of  observation 
array  elements.  When  the  spatial  correlation  ma¬ 
trix  must  be  estimated  from  a  finite  amount  of 
data,  De  Graaf  and  Johnson  (1985)  show  how  the 
maximum-likelihood  estimator  becomes  increas¬ 
ingly  unstable  as  the  data  set  is  reduced. 

IGiapp  and  Carter  (1977)  and  Carter  (1977)  dis¬ 
cussed  the  impact  of  source  motion  on  bias  in 
bearing  estimation  from  a  two-  or  three-element 
array  using  a  time-delay  estimation  procedure. 
Gragg  (1990)  showed  that  there  is  a  time  depen¬ 
dence  in  the  Doppler-shifted  spectra  observed  by 
a  single  receiver  when  an  acoustic  source  has  a 
velocity  vector  that  does  not  pass  through  the 
receiver  location.  He  showed  that  this  time  depen¬ 
dence  makes  the  Doppler  shift  difficult  to  predict. 
Unless  source  motion  effects  can  be  corrected, 
there  will  be  a  phase  bias  in  the  Fourier  estimate  of 
the  time  series.  The  impact  of  the  phase  bias  from 
source  motion  on  wavenumber  estimates  remains 
unquantified. 

In  the  following  sections  we  discuss  the 
beamforming  equations  and  explicitly  define  the 
frequency  domain  estimation  process  used  to  gen- 
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erate  the  cross-power  matrix  in  the  beamforming 
equations.  A  confidence  interval  for  the  variance 
in  the  cross-power  matrix  is  developed.  This  is 
followed  by  a  discussion  of  the  experimental  set¬ 
ting,  data  acquisition  procedure,  and  seismic  ar¬ 
ray  characteristics.  Processing  results  are  then  given 
for  an  impulsive  acoustic  source,  a  pure  seismic 
source,  and  a  moving  M60  tank.  We  obtain  accu¬ 
rate  wavenumber  estimates  from  the  beamformer 
for  all  three  source  types.  Next  we  show  that  small 
variations  in  the  parameters  used  to  estimate  the 
cross-power  matrix  contribute  to  wavenumber 
bias  effects.  These  bias  effects  are  summarized  for 
a  variety  of  block  lengths,  block  overlaps,  and 
window  types. 


2.  BEAMFORMING  PRINCIPLES 

The  discussion  given  here  follows  Capon  (1969) 
and  Johnson  (1982).  First,  a  general  expression  for 
an  F-K  beamformer  power  function  is  given.  The 
standard  Bartlett  (BT)  and  high-resolution  Capon 
maximum-likelihood  (ML)  beamformer  equations 
are  obtained.  Following  the  beamformer  deriva¬ 


tion,  the  estimation  and  conditioning  procedures 
of  the  spatial  correlation  matrix  are  given. 


2.1  Beam-power  response  function 

Assume  an  array  of  M  sensors  have  vector  loca¬ 
tions  Rm  in  the  x-y  plane 


f 

Rm 

=J*m+ym 

V 

> 

as  shown  in  Figure  2.  Each  sensor  is  a  spatial 
wavefield  sample  point  that  records  over  time. 
Such  a  collection  of  sensor  data  is  termed  the  signal 
vector.  A  wavefield  containing  asingle  plane  wave 
source  sampled  at  locations  Rm  may  be  repre¬ 
sented  by  a  signal  vector  with  components: 


*m(0  ~  S 


Rm*  uv 

f  + - -l+nm(f) 


(1) 


where  k0  is  the  wavenumber 


k  o 


=  2tt//c 


/«  k« 


is  the  unit  vector  direction  of  the  source,  nm  is  the 
background  noise  field,  and  c  is  the  scalar  wave 
propagation  velocity.  In  the  frequency  domain, 
eq  1  is 


Figure  2.  A  generic  array  geometry  composed  of  M  sensors  located 
at  vector  positions  Rm.  The  pseudo-plane  wave  source  has  a 
wavenumber  vector  k0 .  The  observation  wavenumber  is  k  .  Dots 
(•)  are  sensor  locations. 
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(8) 


Xm{f)  =  S(/)e'^o  *  Rm  +  Nm(f),  (2) 

where  S(f)  is  the  ideal  Fourier  transform  of  the 
signal  without  the  noise. 

In  matrix  notation,  the  signal  vector  is 


*oCO  \ 

*\(f) 

*2<f) 


x(f}=\  • 


\XM-m 


(3) 


Given  M  spatial  samples  of  the  wave  field,  a 
beam  response  function  [Y(/,  k)]  is  formed  by 
summing  over  each  component  of  the  signal  vec¬ 
tor  and  applying  a  phase  shift  that  is  dependent 
upon  an  observation  wavenumber  ( k) : 

,  ^  M-l 

Y[f,k)=  I  flme-ik  *  Xm  (/) .  (4) 

m  =  0 

Using  eq  2  and  4, 


r(/,fc)=  I  am(f)[s(f)ei*o-Z)-ln  +  Nn 

m  =  0 


ml 

(5) 


where  the  am(f)  constants  are  spatial  filter 
weights  to  be  defined  through  the  application  of  a 
specific  beamforming  strategy,  and  Nm  is  the 

noise  field  modified  by  a  phase  shift  of  elk  *  Rm . 
Equation  5  is  equivalent  to  an  estimation  of  the 
complex  wavenumber  spectrum  of  the  sampled 
wavefield.  In  matrix  notation  the  beam  response  is 

f  \ 


=  A*TX 


(6) 


V  J 


where  *T  means  complex  conjugate  transpose, 
and 


A  (f)  = 


I  flo(f)e-ik*^  \ 

a2{f)e-^ '  ^2 


(7) 


A (f)  is  referred  to  as  the  steering  vector.  The  beam 
power  at  a  given  frequency  and  observation 
wavenumber  is 


P(f,  k)  =  A‘T  XX*T  A 
or 

Pfc  fc)  =  A*TR(/)A,  (9) 

where  R (f)  is  an  M  x  M  matrix  composed  of  the 
product  of  the  signal  vector  and  its  complex  conju¬ 
gate  transpose.  The  matrix  R  is  commonly  referred 
to  as  the  spatial  correlation  matrix.  A  key  result  of 
this  report  will  be  to  show  that  subtle  changes  in 
the  parameters  used  to  estimate  the  correlation 
matrix  affect  the  bias  in  the  beam  power  function. 

Equation  9  is  the  foundation  of  the  wavenum¬ 
ber  estimation  procedure  used  in  this  paper.  Re¬ 
ferring  to  eq  5,  we  see  that  the  maximum  of  the 
beam  power  function  occurs  when  Tcq-T<:=Q. 
When  performing  a  narrowband  F-K  estimation, 
eq  9  is  evaluated  at  a  single  frequency  over  a  reg¬ 
ularly  spaced  grid  of  observation  wave-numbers. 

A  wideband  beam  power  function  can  be  pro¬ 
duced  by  integrating  eq  9  over  a  band  of  frequen¬ 
cies: 


Recently  Nawab  et  al.  (1985)  formulated  more 
elegant  methods  of  computing  a  broadband 
beamformer. 

2.2.  Bartlett  and  Capon  maximum-likelihood 
beamformers 

The  BT  and  ML  methods  are  special  cases  of  eq 
9  and  are  a  consequence  of  defining  specific  values 
for  the  am(f)  weighting  functions  (eq  5).  The  BT 
method  is  the  conventional  method  of  estimating 
an  F-K  spectrum.  The  ML  method  attempts  to 
minimize  energy  leakage  into  the  beam  response 
from  regions  outside  the  observation  wavenumber. 
This  energy  leakage  problem  is  the  spatial  equiva¬ 
lent  of  frequency  domain  bias.  Thus,  am(f)  is  a 
component  of  a  spatial  window  function  analo¬ 
gous  to  time  series  windowing. 

If  all  am(f)  =  1  then  eq  9  yields  the  standard  BT 
beamformer.  Under  this  condition,  eq  8  becomes 


/c-it-So  \ 

!  — *  -* 

e-i  k  •  Ri 

A'  = 

e-ik  •  R2 

• 

• 

1 

\e-i  k  •  rM-1  / 

The  power  function  is 
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PBT(/i)  =  A,TRA',  (12) 

where  R  is  the  estimated  spatial  correlation  ma¬ 
trix. 

In  the  ML  beamformer,  am(f)  are  defined  by  the 
properties  of  the  spatial  correlation  matrix  so  that 
spatial  energy  leakage  is  minimized.  This  is  done 
by  minimizing  the  weighted  array  power  function 
under  the  constraint  that  a  fixed  gain  be  main¬ 
tained  at  the  observation  wavenumber.  Johnson 
(1982)  sets  up  the  problem  as 

=  JL(wTrw)  =  0, 

“VCAr  '  ' 


subject  to  W*T  A'  =  1  .  (13) 

A'  is  given  by  eq  11,  and  W  is  defined  as  the 
weighted  steering  vector  whose  elements  are 

Wm  —  flmA'm  #  (14) 


The  solution  to  eq  13  is  (Johnson  1982,  Capon  1969) 


W  = 


RlA' 
A'*t  R-1A' 


(15) 


source  wavefield  the  bias  is  substantially  smaller 
for  ML  compared  with  the  BT  method  (De  Graaf 
and  Johnson  1985). 


2.3.  Estimation  and  conditioning  of  the 
correlation  matrix 

The  spatial  correlation  matrix  (R)  is  constructed 
from  the  frequency  domain  representation  of  the 
signal  vector  x.  The  properties  of  the  spatial  corre¬ 
lation  matrix  determine  the  quality  of  the 
beampower  function.  To  estimate  R  accurately, 
we  use  a  tapered  overlapped  block-averaging  FFT 
(OBAFFT).  This  is  a  widely  utilized  procedure  for 
estimating  signals  with  intricate  spectra  (Welch 
1967,  Capon  et  al.  1967).  The  OBAFFT  method 
allows  one  to  make  trade-offs  between  the  bias 
and  variance  in  the  estimate  of  R. 

It  is  important  to  guarantee  the  nonsingularity 
of  the  spatial  correlation  matrix  when  applying 
the  ML  processing  procedure.  This  can  be  done  by 
the  block  averaging  process  or  it  can  be  forced  by 
applying  a  small  amount  of  additive  noise  to  the 
correlation  matrix  (Capon  1969).  A  normalization 
procedure  that  produces  ones  along  the  diagonal 
of  the  spatial  correlation  matrix  is  also  useful  in 
making  qualitative  assessments  of  the  accuracy  of 
the  wavenumber  spectrum. 


Note  that  eq  15  is  a  function  of  the  spatial  correla¬ 
tion  matrix.  The  effect  of  designing  the  weighted 
steering  vector  using  the  observed  signal  proper¬ 
ties  is  that  we  produce  a  wavenumber  window 
function  that  is  adaptive  to  the  signal-in-noise 
field  and  the  specific  array  geometry.  The  result¬ 
ant  beampower  function  is 

PmlW  =  W*Trw.  (16) 

Substituting  eq  15  into  16  and  using  the  estimation 
for  the  cross-power  matrix,  we  obtain 

PmlM=  - ,  (17) 

A'*TR“V 

where  A  indicates  an  estimate  based  on  available 
data. 

This  is  the  working  form  of  the  ML  beamformer 
used  in  this  paper.  It  should  also  be  noted  that  the 
ML  beamformer  is  also  termed  a  minimum  vari¬ 
ance  method,  and  in  the  case  of  only  one  coherent 
signal  it  is  also  spatially  unbiased  (Capon  1969),  in 
the  sense  that  the  bias  tends  toward  zero  as  the 
number  of  sensors  becomes  large.  In  a  multiple- 


2.3.1.  Tapered  overlapped  block-averaged  FFT 

The  mathematical  procedure  used  to  form  the 
OBAFFT  is  as  follows.  The  time-domain  signal 
vector  obtained  from  a  passive  array  with  M  sen¬ 
sors  and  components  such  as  that  of  eq  1  has  the 
form 


*o(0  \ 
*1  (0 
*2  (0 


\*M-1  (0/ 


Let  x  have  N  samples  at  intervals  of  At.  Now 
define  a  block  of  this  sequence  with  a  length  L, 
such  that  L<N ;  further,  let  each  block  of  L  samples 
be  overlapped  by  a  percentage  p.  The  number  of 
blocks  will  be 

nblks  =  1  +  [— — — ] ,  (19) 

\nshftl 

where  nshft  -  INT[L  x  (l  -  p)]  and  INT  is  the  near¬ 
est  integer  operator.  Thus,  nblks  gives  the  total 
number  of  overlapped  segments  (or  blocks)  that 
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(21) 


are  overlapped  by  an  integer  number  of  points. 
The  Fourier  transform  of  the/th  block  from  chan¬ 
nel  m  will  be 

if )  =  7  X  atxm  (*  +  ]•  nshft)e~l2nlm) 

1  i=o  (20) 

where  the  a ^  are  time-domain  window  taper 
weights.  The  window  tapers  decrease  the  bias  in 
the  estimate  by  ensuring  that  the  ends  of  each 
segment  of  data  are  gradually  tapered  to  zero. 

A  variety  of  window  taper  functions  may  be 
applied,  each  having  a  different  impact  on  the 
properties  of  the  resulting  spectral  estimate.  Har¬ 
ris  (1978)  presents  an  excellent  overview  on  win¬ 
dow  tapering  strategies  for  harmonic  analysis 
problems.  Figure  3  graphically  illustrates  the 
overlapped  (OBAFFT)  segmentation. 

2.3.2.  Formation  of  spatial  correlation  matrix 

The  cross-power  function  between  all  pairs  of 
sensors  in  the  array  is  calculated  for  each  block  of 
data  in  the  OBAFFT  method.  The  final  estimate  of 
the  spatial  correlation  matrix  is  formed  by  averag¬ 
ing  over  the  block  cross-power  matrices. 

Using  eq  20,  elements  of  the  cross-power  matrix 
for  the  ;th  block  are 


<n(f)  =  ^ln(/)*<i(/), 

where  m  and  n  are  sensor  indices.  The  final  esti¬ 
mate  of  the  spatial  correlation  matrix  is  formed  by 
averaging  over  the  block  index, 

Vn(f)  =  —  £<n  •  (22) 

nblks  ,=1 

The  formation  of  the  cross-power  function  be¬ 
tween  two  vectors  produces  a  vector  that  has  a 
phase  angle  that  is  the  difference  in  phase  between 
the  two  constituent  vectors.  For  waveforms  pro¬ 
duced  by  a  stationary  source  and  recorded  by  chan¬ 
nel  m,  the  Fourier  transform  at  any  given  fre¬ 
quency  will  be  a  randomly  oriented  vector.  Thus, 
the  OBAFFT  gives,  for  each  block  of  data  in  chan¬ 
nel  m,  a  series  of  vectors  with  random  phase  angles. 
In  an  array  of  sensors  recording  real  data,  the  quan¬ 
tity  to  be  measured  is  the  phase  difference  between 
channels.  The  final  estimate  of  the  spatial  correla¬ 
tion  matrix  given  by  eq  22  is  therefore  an  average 
of  several  phase  angle  difference  matrices,  which 
becomes  less  variable  as  the  number  of  blocks 
increases.  This  assumes  that  the  phase  difference 
between  channels  is  not  changing  with  time. 


Time  (s) 
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Figure  3.  Segmentation  and  overlapping  of  a  time-domain  component  of  the 
signal  vector.  Each  tapered  block  is  transformed  into  the  frequency  domain  via 
an  FFT. 


6 


X0(/)XoV) 

Xi(/)X0V) 

X2(/)X0V) 

x„  wl(f) 

Wi(/) 

X2(f)X*(f) 

x„  wKv) 
Xj(/)X2V) 
X2(/)X2V) 

•  x0c 

•  XiWu-iif) 

■  X2(/)xLi  if) 

R(Z)  = 

• 

• 

• 

xM-i  ml(f) 

V 

• 

• 

• 

XM-l(/)XiV) 

• 

• 

• 

XM-i(/)X2V)  • 

• 

• 

• 

) 

(23) 

Figure  4.  Spatial  correlation  matrix  structure. 


2.3.3.  Spatial  correlation  matrix  conditioning 

The  structure  of  the  spatial  correlation  matrix  is 
given  by  eq  23  (Fig.  4).  In  the  programmed  imple¬ 
mentation,  each  element  of  this  matrix  is  normal¬ 
ized  according  to 


/V  A* 


XjCQXjCQ 

|xid|xjd 


(24) 


If  the  signal  vector  is  perfectly  coherent,  then  the 

diagonal  elements  of  R  will  contain  unity  values, 
resulting  in  a  maximum  beam  power  response  of 
1  when  applying  the  BT  method.  This  allows  quali¬ 
tative  judgments  to  be  made  on  the  reliability  of 
the  wavenumber  estimate,  which  is  particularly 
important  when  dealing  with  nonideal  signal  types 
or  inhomogeneous  propagation  environments. 
When  using  narrowband  BT  processing  with  this 
normalization,  a  peak  power  response  <  0.5  is 
often  unreliable. 

The  existence  of  the  R  matrix  is  of  critical 
importance  when  applying  the  ML  method  (see  eq 
18).  Capon  (1969)  notes  that  if  the  number  of  blocks 
(N)  in  the  estimate's  less  than  the  number  of  array 
sensors  (M),  then  R  is  of  order  M  and  rank  N  and 
is  thus  singular.  This  is  frequently  the  case  when 
processing  signals  of  short  duration.  To  guarantee 
the  nonsingularity  of  the  estimated  spatial  corre¬ 
lation  matrix  we  add  a  small  amount  of  incoherent 
noise  (X)  to  the  elements  of  the  spatial  correlation 
matrix.  Capon  (1969)  suggests 

R'  =  (l- A,)R  +M.  (25) 


The  choice  of  X  is  determined  by  trial  and  error  (the 
results  reported  here  often  used  values  <  1C h4).  In 
practice  the  accuracy  of  the  beam  response  is  rela¬ 


tively  insensitive  to  X,  but  it  should  be  kept  as 
small  as  possible  since  it  lowers  the  signal-to-noise 
ratio  (SNR)  of  the  beam  response. 

2.3.4.  Variance  of  the  cross-power  matrix 

The  OB AFFT  method  allows  a  degree  of  control 
over  the  bias/variance  trade-off  in  the  estimate  of 
the  cross-power  matrix.  The  bias  is  mitigated  by 
the  choice  of  window  taper  and  the  number  of 
points  in  each  block,  while  the  variance  is  a  func¬ 
tion  of  the  number  and  degree  of  overlap  between 
blocks.  The  independence  of  each  block  is  reduced 
by  increasing  the  overlap.  The  characteristics  of 
the  window  taper  applied  to  each  block  also  affect 
the  block's  independence  for  a  given  degree  of 
overlap.  This  partial  dependence /independence 
may  be  exploited  to  decrease  the  variance  of  a 
given  spectral  estimate. 

As  an  illustration  of  the  variance  reduction  that 
can  be  obtained  using  this  strategy,  a  signal  (dt  = 
1  / 1023)  was  generated  by  superimposing  two  sine 
waves  with  line  spectra  at  20  and  23  Hz  and 
unitary  peak-to-peak  amplitudes.  The  second  sine 
wave  series  is  time-shifted  (relative  to  the  first 
signal)  by  0.0455  s  (at  20  Hz  in  the  frequency 
domain  this  is  an  approximately  32°  phase  shift). 
A  unique  pseudo-random  Gaussian  "noise"  series 
with  zero  mean  and  a  maximum  peak-to-peak 
amplitude  of  1  is  added  to  each  of  the  two  sine 
wave  signals.  This  yields  an  SNR  of  1  for  each 
spectral  line  in  the  series.  The  OBAFFT  spectral 
estimation  process  was  applied  at  20  Hz  for  each 
sine  wave  plus  noise  signal,  and  the  cross-power 
operation  was  carried  out.  Note  that  the  presence 
of  the  23-Hz  spectral  line  in  the  resulting  series 
allows  simulation  of  phase  and  amplitude  bias 
effects.  The  20-Hz  OBAFFT  estimate  of  the  cross 
power  between  the  two  noisy  signals  was  ob- 


7 


Estimate  of  Normalized  Cross-Power  Amplitude 
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Estimated  Cross-Power  Phase  Values  (degrees) 


Estimate  of  Normalized  Cross-Power  Amplitude  Estimated  Cross-Power  Phase  Values  (degrees) 


Figure  5.  Histograms  of  cross-power  amplitude  and  phase  estimates  at  20  Hz  for  a  synthetic  sine  wave  series 
containing  line  spectra  at  20  and  23  Hz  with  zero  mean  random  noise.  SNR  is  one.  Sample  populations  for  all 
histograms  are  composed  of 5000  estimates.  All  estimates  utilized  a  Blackman  window  taper  to  reduce  bias  effects  from 
23  Hz  spectral  line.  Cross-power  histograms  shown  in  (a)  and  (b)  are  based  on  estimates  from  blocks  of  1024  data 
points.  Each  estimate  in  the  amplitude  and  phase  histograms  shown  in  (c)  and  (d)  are  based  on  six  overlapped  blocks 
of  512  points  each.  The  total  series  length  is  1024  points ,  and  each  block  overlaps  the  preceding  block  by  80%.  Note 
that  the  variance  in  histograms  (c)  and  (d)  is  significantly  smaller  than  the  variance  seen  in  (a)  and  (b). 


tained  for  5000  different  random  noise  cases. 
Using  this  sample  population  of  cross-power  val¬ 
ues,  the  mean  and  variance  of  the  amplitude  and 
phase  population  can  be  estimated. 

Figure  5  shows  cross-power  amplitude  and 
phase  estimate  histograms  using  the  Blackman 
window  function  (given  by  W(l)  -  0.42  -  0.5  cos 

(r^r)  +  cos  (r^r)  /  =  1,  L).  Figure  5a 

and  5b  are  the  amplitude  and  phase  histograms  for 
the  case  of  estimates  based  on  a  single  block  of  1024 
points  tapered  with  a  Blackman  window  function. 


There  are  two  degrees  of  freedom  in  this  cross¬ 
power  estimate.  Based  on  a  sample  population  of 
5000  estimates,  the  amplitude  mean  is  1.014  with  a 
variance  of  0.0829.  The  phase  mean  is  39.36  °  with 
a  variance  of  4.731°2.  In  Figure  5c  and  5d  the 
amplitude  and  phase  histograms  are  calculated 
but  the  series  is  estimated  using  blocks  overlapped 
by  80%.  Thus,  six  partially  independent  blocks  are 
used  in  each  amplitude  and  phase  estimate.  Again 
using  a  population  of  5000  estimates,  the  histo¬ 
gram  shown  in  Figure  5c  gives  an  amplitude  mean 
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of  1.013  with  a  variance  of  0.05711.  In  Figure  5d  the 
phase  estimate  population  yields  a  mean  of  39.75° 
with  a  variance  of  3.249°2.  The  critical  point  to  note 
is  the  decrease  in  the  variance  for  the  case  of  the 
overlapped  six  block  estimates  as  compared  with 
the  variance  calculated  for  the  single  block  esti¬ 
mates.  Thus,  it  is  statistically  advantageous  to  use 
OBAFFT  to  increase  the  stability  of  and  decrease 
the  bias  in  spectral  estimation. 

2.4.  Verification  of  signal-processing  algorithm 
The  Fortran  code  implementing  the  beamform¬ 
ing  procedures  given  above  is  discussed  in  detail 
by  Moran  (1991).  The  algorithm  was  extensively 
tested  during  development  and  continues  to  be 
refined.  Two  tests  were  used  to  validate  the  proce¬ 
dure.  The  first  test  was  numerical,  using  a  wide 
variety  of  synthetic  signals  under  various  signal- 
to-noise  ratios  (SNRs)  and  with  a  variety  of  fre¬ 
quency  bands  and  array  geometries.  The  second 
and  most  critical  test  came  from  results  obtained 
from  field  data  collected  under  defined  experi¬ 
mental  conditions.  These  field  data  tests  are  based 
on  the  impulsive  sources  discussed  in  section  5.2. 


Figure  6.  Location  map  for  the  seismic  experiments , 
showing  vehicle  test  track ,  location  of  contractor's  sen¬ 
sors, and  location  ofCRREL  seismic  arrays.  Also  shown 
is  the  location  of  the  snow  pits ,  met  station,  frost  tubes , 
and  borehole. 


3.  FIELD  DATA  ACQUISITION  AND 
EXPERIMENTAL  CONDITIONS 

To  test  the  efficacy  of  beamforming,  experimen¬ 
tal  measurements  of  various  sources  were  ob¬ 
tained.  The  environmental  characteristics  of  the 
test  site  and  the  methods  used  to  make  these 
recordings  are  discussed  in  this  section. 

3.1.  Geological  setting 

The  field  tests  were  performed  in  February  1988 
at  Range  37  at  the  Michigan  National  Guard's 
Camp  Grayling  (Fig.  6).  The  soil  in  the  Grayling 
locale  consists  of  glacial  outwash  sand  and  gravel 
(Farrland  and  Bell  1982).  Logs  from  four  wells  in 
the  vicinity  show  the  thickness  of  these  deposits  at 
116, 137, 175,  and  251  m,  with  the  bedrock  below 
being  Mississippian  in  age  (Lilienthal  1978).  A 
local  driller  estimated  the  sand  thickness  in  the 
area  to  be  around  180  m  based  on  his  experience. 
A  borehole  30  m  deep  was  drilled  at  the  site;  the 
soils  logged  during  the  drilling  are  given  in  Table 

1.  The  characteristics  of  eight  soil  samples  that 
were  collected  off  the  drill  auger  are  listed  in  Table 

2.  These  observations  show  that  the  site  was  un¬ 
derlain  by  medium  sand,  with  some  thin  layers  of 
fine  gravel  and  clay  at  depths  around  4  and  24  m. 
The  water  table  was  5.8  m  below  the  surface. 

3.2.  Shallow  seismic  velocity  structure 

Standard  compressional  (P)  and  shear  (S)  wave 

refraction  methods  were  applied  at  the  site  to 
determine  the  seismic  velocity  structure.  A  linear 
array  of  24  geophones  spaced  1  m  apart  was  ar¬ 
ranged  perpendicular  to  the  test  track.  Hammer 
blows  were  then  recorded  at  intervals  off  each  end 
of  the  array.  The  resulting  distance  vs.  travel  time 
curves  were  analyzed  using  the  intercept  time 
method  to  determine  the  seismic  velocity  as  a 
function  of  depth.  For  the  P  waves,  vertical  sledge¬ 
hammer  blows  on  a  metal  plate  at  the  ground 
surface  were  used  as  the  source,  and  Mark  Prod¬ 
ucts  Model  L15-B  vertical-component  geophones 
with  resonant  frequencies  of  4.5  Hz  were  used  to 
record  the  resulting  ground  motion.  For  the  S 
waves,  a  0.2  x  0.2  m  wooden  plank  was  aligned 
perpendicular  to  the  array  and  clamped  to  the 
ground  by  the  front  wheels  of  a  pickup  truck  (Fig. 
7).  Horizontal  sledgehammer  blows  on  each  end 
of  the  plank,  which  excited  SH  (shear  horizontal) 
waves  of  opposite  polarity,  were  then  recorded  at 
each  source  location  using  horizontal-component 
geophones. 

One  difficulty  encountered  was  caused  by  the 
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Table  1.  Borehole  soil  log.* 


Depth  to  Thickness 

bottom  of  layer  of  layer 

Layer  (m)  (ft)  (m)  (ft)  Soil  type 


1 

3.4 

11 

3.4 

11 

Medium  sand 

2 

4.9 

16 

1.5 

5 

Fine  gravel 

3 

5.8 

19 

0.9 

3 

Medium  sand 

4 

7.6 

25 

1.8 

6 

Medium  sand  and  water 

5 

7.9 

26 

0.3 

1 

Clayish  fine  sand 

6 

23.2 

76 

15.2 

50 

Medium  sand  and  water 

7 

23.5 

'  77 

0.3 

1 

Blue  clay 

8 

23.8 

78 

0.3 

1 

Medium  and  sand  water 

9 

24.4 

80 

0.6 

2 

Blue  clay 

10 

24.7 

81 

0.3 

1 

Brown  clay 

11 

27.4 

89 

2.7 

9 

Medium  sand  and  water 

12 

29.3 

96 

1.8 

6 

Clean  sand  and  water 

13 

30.5 

100 

1.2 

4 

Medium  sand  and  water  with 

clay  trace  in  water 


*From  Hal  Carlson,  Jim's  Well  Drilling,  Grayling,  Michigan. 


presence  of  the  frost  layer.  For  recordings  made 
close  to  the  hammer  source,  the  first  arrival  was  a 
wave  traveling  directly  through  the  frost  layer,  at 
a  speed  of  about  590  m  s-1 .  Fortunately,  this  arrival 
was  highly  attenuated  and  died  out  within  a  few 
meters.  We  were  able  to  pick  the  travel  times  of  the 
later-arriving  refracted  waves  after  this  first  ar¬ 
rival  and  use  them  in  the  analysis.  The  P-wave 
velocity  of  the  thin  snow  layer  was  measured  at 
about  100  m  s-1  by  a  short  line  of  geophones  placed 
above  and  below  the  snow  cover  and  spaced  0.25 


m  apart.  The  measured  first  arrival  travel  times  for 
the  P-  and  S-wave  refraction  experiments  along 
with  least-squares  line  fits  are  shown  in  Figures  8 
and  9. 

Analysis  of  the  travel  times  indicates  a  surface 
layer  4.8  m  thick  with  a  P-wave  velocity  of  290 
m  s*-1,  underlain  by  a  layer  with  a  velocity  of  1660 
m  s^1  that  is  identified  as  the  water  table.  The  shear 
wave  data  reveal  a  surface  layer  with  a  velocity  of 
150  m  s”1, 12  to  15  m  thick,  sloping  downward  at  an 
angle  of  about  9  degrees  to  the  north,  above  a  layer 


Table  2.  Physical  properties  of  borehole  soil  samples. 


Percentage  of  soil  type 


Moisture 


Sample 

no. 

Depth 
(m)  (ft) 

Description 

Density 
(kg  m3) 

content 

(%) 

Gravel 

fine 

Sand 

coarse 

Sand 

medium 

Sand 

fine 

Clay 

1 

1.5 

5 

Sand 

1.38 

3.7 

7.6 

1.8 

34.5 

63.5 

0.3 

2 

3.0 

10 

Sand 

1.45 

3.4 

4.6 

2.0 

33.0 

63.9 

1.0 

3 

4.0 

13 

Sand  and  gravel 

1.56 

2.8 

18.2 

10.2 

35.7 

35.6 

1.1 

4 

6.1 

20 

Sand  (below 
water  table) 

1.45 

3.4 

0.4 

0.3 

12.9 

86.7 

0.4 

5 

13.7 

45 

Sand 

1.45 

19.2 

0.0 

0.2 

10.0 

88.1 

1.7 

6 

22.9 

75 

Blue  clay 

— 

23.2 

0.0 

0.1 

7.1 

53.4 

39.4 

7 

24.4 

80 

Brown  clay 

— 

6.8 

0.0 

0.1 

4.4 

35.5 

60.0 

8 

30.5 

100 

Sand  with  clay 

1.73 

19.4 

0.0 

0.0 

15.0 

75.7 

9.4 

The  percentage  of  gravel,  sand,  and  clay  in  each  sample  was  determined  from  the  grain  size  measurements.  According  to  the  unified 
soil  classification  scheme,  soils  with  grains  smaller  than  9.5  mm  but  larger  than  2  mm  are  classified  as  fine  gravel.  Sand  has  grain  sizes 
less  than  2  mm  and  greater  than  0.076  mm  and  is  further  divided  into  coarse,  medium,  and  fine  at  the  0.59  mm  and  0.25  mm  sizes. 
Grains  less  than  0.076  mm  can  be  either  silt  or  clay,  but  are  assumed  here  to  be  clay  based  on  field  observations  of  the  samples. 
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Figure  7.  Generating  shear  waves  for  the  site  characterization  measurements.  The  cable  leading 
from  the  sledgehammer  allows  the  impact  time  to  be  recorded. 


Figure  8.  Compressional  (P)  wave  refraction  travel 
time  vs.  distance.  The  straight  lines  are  least-squares 
fits  to  the  travel  time  data.  Shown  are  arrival  times  for 
waves  from  hammer  sources  located  north  and  south  of 
the  linear  array  of  geophones. 


Figure  9.  Shear  (S)  wave  refraction  travel  time  vs. 
distance.  The  straight  lines  are  least-squares  fits  to  the 
travel  time  data.  Shown  are  arrival  times  for  waves  from 
hammer  sources  located  north  and  south  of  the  linear 
array  of  geophones. 
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Table  3.  Analysis  of  seismic  refraction  travel  time. 


Apparent  velocity  Intercept  time  Depth  Dip 

(m  s-i)  (ms)  Velocity  (m)  angle 


Layer 

South 

North 

South 

North 

(m  s~J) 

South 

North 

(degrees) 

Compressional  (P)  waves 

1 

276 

301 

0.0 

0.0 

288 

4.78 

4.80 

2 

1584 

1749 

32.6 

32.8 

1662 

-0.5 

Shear  (S)  waves 

1 

167 

137 

0.0 

0.0 

152 

12.2 

15.7 

2 

417 

386 

148.3 

191.2 

401 

with  a  speed  of  400  m  s-1.  Details  of  the  intercept 
time  analysis  are  given  in  Table  3.  The  seismic 
water  table  depth  is  about  1  m  less  than  the  well  log 
depth,  while  the  shear  wave  velocity  change  does 
not  correspond  to  any  obvious  drill  log  layer.  The 
vertical  compressional  and  shear  wave  velocities 
of  the  shallow  soil  were  also  measured  by  clamp¬ 
ing  a  three-component  geophone  at  various  depths 
in  a  cased,  30-m-deep  borehole  and  recording  the 
travel  times  of  P  and  S  waves  from  surface  hits  3  m 
away  from  the  top  of  the  hole  to  the  downhole 
geophone.  The  seismograms  from  the  two 


Figure  10.  Compressional  (P)  wave  downhole  travel 
time  vs.  distance.  Straight-line  ray  paths  were  assumed 
to  calculate  the  slant  distance  from  the  measured  sensor 
depth  and  horizontal  source  offset.  These  lines  yielded 
values  of 390  m  s~2  down  to  3  m,  underlain  by  a  velocity 
of  1370  m  s-J. 


downhole  horizontal  components  of  unknown 
orientation  were  rotated  to  the  radial  and  trans¬ 
verse  orientations  before  the  first  arrivals  were 
picked.  Figures  10  through  12  show  the  resulting 
travel  time  curves  for  the  P  and  S  waves,  with 
least-squares  line  fits  used  to  determine  the  veloci¬ 
ties.  These  experiments  indicate  compressional 
wave  velocities  of  390  m  s-1  at  the  surface,  with  a 
1370-m  s-1  layer  below,  starting  at  a  depth  of 
around  3  m.  The  measured  shear  wave  velocity 
was  180  m  s_1  at  the  surface  and  340  to  380  m  s-1  at 
a  depth  of  12  to  16  m.  These  results  agree  with  the 


Figure  11.  Transverse  shear  (S)  wave  downhole  travel 
time  analysis.  The  least-squares  lines  give  180  m  s -1 
down  to  12  m  and  340  m  s~ 1  below  that  depth. 
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Figure  12.  Radial  shear  (S)  wave  downhole  travel  time 
analysis.  A  velocity  of  180  m  s-1  down  tol6mf  with  380 
m  s-1  below  that ,  is  indicated  by  the  least-squares  line 
At. 


refraction  results  discussed  above.  The  average 
velocities  over  the  upper  30  m  were  1090  m  s_1  for 
the  P  waves  and  250  m  s_1  for  the  two  S  waves. 

3.3.  Environmental  characterization 

A  number  of  meteorological  parameters  were 
recorded  continuously  at  the  site  of  the  seismic 
and  acoustic  experiments.  A  Campbell  Scientific 
Inc.  (CSI)  Model  21X  Micrologger  was  used  to 
collect  the  data  (Fig.  13),  which  were  recorded  on 
a  cassette  magnetic  tape  and  written  on  a  paper 
tape.  CSI  Model  107  thermistors  recorded  the  tem¬ 
peratures  in  the  ground,  in  the  snow,  and  at  heights 
of  0.01, 0.1, 0.3, 1,  and  1.4  m  above  the  snow.  These 
readings  have  an  accuracy  of  ±0.2°C.  A  CSI  Model 
207  sensor  provided  the  temperature  and  relative 
humidity  at  a  height  of  2  m.  The  wind  speed  and 
direction  at  a  height  of  3  m  was  recorded  by  Met 
One  Model  014A  and  024A  sensors  with  an  accu¬ 
racy  of  ±2%  and  ±5%,  and  a  YSI  Model  2014  barom¬ 
eter  provided  the  barometric  pressure  to  within  1 
mb.  Readings  from  the  12  sensors  were  taken 
every  minute,  and  the  average  and  standard 


Figure  13.  The  instrumentation  used  to  record  meteorological  pa¬ 
rameters.  The  Campbell  micrologger  is  housed  in  the  box  at  the  center 
of  the  tripod.  At  the  top  are  the  wind  speed  and  direction  sensors.  The 
radiation  shield  contains  a  thermistor  and  the  relative  humidity 
sensor.  Additional  thermistors  and  a  barometer  are  mounted  on  the 
wooden  post.  A  lead-acid  battery  below  the  tripod  provided  the 
power. 
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Table  4.  Meteorological  data. 

Temperatures  (°C) 


Date 

Time 

Gnd 

Snow 

0.02  m 

0.2  m 

0.3  m 

2  m 

1.4  m 

2  m 

Pressure 

(mb) 

RH 

(%) 

WS 

(ms-D 

WD 

(deg) 

16  Feb 

1630 

-3.4 

-3.4 

-3.9 

-4.2 

-4.4 

-4.2 

-3.5 

-4.5 

963.6 

94 

3.3 

201 

16  Feb 

1700 

-3.3 

-3.4 

-4.0 

-4.2 

-4.4 

-4.2 

-3.5 

-4.4 

963.8 

93 

4.0 

209 

16  Feb 

1730 

-3.3 

-3.4 

-5.2 

-5.0 

-5.2 

-5.0 

-3.7 

-4.9 

963.9 

92 

3.9 

207 

16  Feb 

1800 

-3.3 

-3.4 

-5.4 

-5.1 

-5.3 

-5.2 

-4.1 

-5.1 

963.9 

92 

3.0 

203 

16  Feb 

1830 

-3.3 

-3.4 

-5.5 

-5.1 

-5.4 

-5.3 

-4.4 

-5.2 

964.0 

92 

2.7 

222 

16  Feb 

1900 

-3.3 

-3.4 

-5.5 

-5.2 

-5.5 

-5.4 

-4.6 

-5.3 

964.0 

93 

3.0 

195 

17  Feb 

1100 

-2.7 

-2.7 

-1.8 

-1.2 

-1.3 

-0.8 

0.0 

-1.5 

971.1 

80 

3.1 

267 

17  Feb 

1130 

-2.7 

-2.7 

-1.9 

-1.1 

-1.1 

-0.8 

0.3 

-1.4 

971.5 

80 

3.5 

311 

17  Feb 

1200 

-2.7 

-2.7 

-2.2 

-1.7 

-1.8 

-1.5 

-0.2 

-1.7 

971.5 

81 

5.1 

306 

17  Feb 

1230 

-2.6 

-2.6 

-1.2 

-0.7 

-0.7 

-0.2 

0.1 

-1.4 

971.8 

76 

5.3 

286 

17  Feb 

1300 

-2.6 

-2.5 

-2.2 

-1.5 

-1.7 

-1.2 

0.4 

-1.4 

971.8 

76 

4.5 

302 

17  Feb 

1330 

-2.5 

-2.5 

-2.5 

-1.6 

-1.6 

-1.2 

0.7 

-1.6 

971.7 

74 

6.4 

256 

17  Feb 

1400 

-2.5 

-2.4 

-1.2 

-0.3 

-0.1 

0.3 

0.9 

-1.2 

971.8 

74 

4.7 

272 

deviation  of  the  values  were  computed  and  re¬ 
corded  every  30  minutes.  The  meteorological  data 
recorded  on  16  and  17  February  1988  during  the 
seismic  and  acoustic  array  experiments  discussed 
in  this  report  are  tabulated  in  Table  4.  Two  frost 
tubes  were  installed  at  the  site  in  January:  one  was 
located  5  m  north  of  the  met  station  in  a  topo¬ 
graphic  high,  and  the  second  was  located  10  m 
south  of  the  met  station  in  a  topographic  low.  The 
elevation  difference  between  the  two  was  0.95  m. 
The  snow  depth  at  the  first  frost  tube  was  0.14  m, 
and  frost  depth  exceeded  the  depth  of  the  tube 
(0.775  m)  for  all  of  the  February  readings.  On  16 
February,  the  snow  depth  at  the  second  tube  was 
0.18  m  and  the  frost  depth  was  0.60  m.  Snow 
depths  at  the  site  during  these  experiments  ranged 
from  0.14  to  0.20  m  in  most  locations.  The  snow 
had  been  on  the  ground  for  some  period  of  time.  A 


snow  pit  examined  on  18  February  revealed  the 
presence  of  a  thin  wind  crust  at  the  top  underlain 
by  rounded  snow  crystals  with  a  layer  of  depth 
hoar  at  the  bottom  (Table  5). 

3.4.  Array  description  and  recording  procedures 

All  the  waveform  data  were  digitally  recorded 
by  a  Geosource  DSS-10A  system  on  9-track  mag¬ 
netic  tape  at  1600  bpi  in  the  SEG  B  format.  The 
recording  system  has  a  dynamic  range  of  90  dB 
and  a  bandwidth  from  3  to  500  Hz  at  a  sampling 
rate  of  2  kHz.  Mark  Products  L-15B  vertical-  and 
horizontal-component  geophones  with  a  resonant 
frequency  of  4.5  Hz  and  a  sensitivity  of  32  V  m 
s-1  were  used  to  detect  the  ground  motion.  Globe 
100C  low-frequency  microphones  with  a  sen¬ 
sitivity  of  2  V  Pa”1  were  also  used  as  sensors.  The 
geophones  and  microphones  were  placed  in  an  18- 


Table  5.  Snow  characterization,  pit  2,  at  1020  on  18  February  1988. 


Layer 

Depth 

(mm) 

Density 
(kg  m~3) 

Tempera¬ 

ture 

(°C) 

Hardness 

index 

Crystal 

size 

(mm) 

Crystal  type 

1 

0 

-3.5 

Wind  crust  (irregular,  broken 

crystals),  planar  dendrites 

1 

120 

-3.5 

60 

0.5 

Rounded,  often  branched 

35 

150 

-5.0 

5 

0.5 

Rounded,  often  branched 

85 

210 

-6.5 

10 

0.5 

Rounded,  often  elongated 

2 

115 

-5.0 

0 

2.5 

Depth,  hoar,  angular. 

stepped  grains 

140 

160 

-4.5 

2.5 

Depth,  hoar,  angular, 

stepped  grains 

165 

-4.5 

2.5 

Depth  hoar,  angular. 

stepped  grains 
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Figure  14.  Grayling  sensor  positions  and  road ,  drawn  to 
scale.  Relative  sensor  coordinates  and  associated  channel 
numbers  are  given  in  Table  6.  North  is  referenced  to  the 
magnetic  pole.  Dots  (•)  are  sensor  locations. 


element  crossed  array  with  one  of  the  legs  perpen¬ 
dicular  to  the  test  track.  Figure  14  shows  the  sensor 
positions,  and  Table  6  gives  the  relative  element 
positions.  The  array  was  calibrated  by  vertical 
hammer  blows  and  blank  pistol  shots  from  five 
locations  along  the  test  track:  at  the  closest  point  of 
approach  of  the  vehicle  (CPA),  and  at  50  and  100  m 
on  either  side  of  the  CPA.  Shear  waves  from  the 
CPA  were  also  recorded,  as  were  pistol  shots  50  m 
from  the  array  center  at  the  remaining  compass 
points. 


Table  6.  Seismic  array  element  position. 


Element 

EtoW(m) 

N  to  S  (m) 

1 

0.0 

-2.0 

2 

0.0 

0.0 

3 

30.0 

6.0 

4 

0.0 

12.0 

5 

0.0 

16.0 

6 

1.73 

-1.0 

7 

-5.0 

3.0 

8 

-10.4 

6.0 

9 

-13.8 

8.0 

3.5.  Characterization  of  the  Grayling  array 
An  important  factor  in  qualitatively  assessing 
the  reliability  of  an  estimated  wavenumber  spec¬ 
trum  based  on  the  BT  or  ML  methods  is  the  pres¬ 
ence  of  the  array  factor  [F{k)]  in  the  estimate.  The 
array  factor  is  the  zero-phase  impulse  response  for 
a  given  array  geometry.  The  array  impulse  re¬ 
sponse  provides  important  information  on  the 
direction-dependent  resolution  capability  of  the 
array  as  well  as  an  estimate  of  the  beam  width  for 
an  ideal  plane  wave.  It  is  a  relatively  simple  matter 
to  obtain  such  a  response.  The  most  straightforward 
method  is  given  by 

M-l 

F(k)=~  X  e-ik'Rm,  (31) 

M  m  =  o 

which  is  simply  a  Fouriertransform  of  the  array 
element  location  vectors  Rm.  The  procedure  used 
in  this  analysis  is  given  by  Kelley  (1967)  and  differs 
from  eq  13  only  in  that  it  allows  for  the  analysis  of 
a  wideband  array  response  based  on  a  Gaussian 
distribution  of  zero-phase  energy. 
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East-West  Component  of  Wavenumber  (nr1) 


Linear  grey  scale:  White  Regions  =  -3  dB 
Black  Regions  =  -42  dB 

Figure  15.  Grayling  array  impulse  response.  Plot  ranges  over  kx  and  ky 
Nyquist  limits.  Spectral  resolution  at  -3  dB  is  0.19  ra~T 


The  impulse  response  for  the  1989  Grayling 
array  is  shown  in  Figure  15.  The  plot  ranges  to  the 
Nyquist  wavenumber  limits  along  the  x  and  y 
axes.  The  array  factor  is  the  x  described  by  the  -12- 
dB  high-amplitude  ridge  lines  radiating  from  the 
center  of  the  plot.  The  inner  contour  at  -3  dB 
defines  the  spectral  resolution  for  this  array,  which 
is  0.19  m-1.  The  high-amplitude  lobes  seen  intrud¬ 
ing  from  the  periphery  of  the  plot  toward  the 
center  are  grating  lobes,  which  are  a  consequence 
of  the  periodic  nature  of  the  spectral  estimation 
procedure. 

The  Grayling  array  geometry  is  not  optimal 
with  respect  to  multiple  source  resolution  or  pro¬ 
viding  equal  resolution  capabilities  from  all  bear¬ 
ing  directions.  Signals  arriving  from  directions  in 
the  vicinity  of  0°  or  180°  and  55°  or  235°  will  be 
poorly  resolved.  Maximum  source  resolution  lies 
at  directions  perpendicular  to  these  values.  The 
element  spacing  of  the  Grayling  array  has  been 
designed  to  allow  coverage  of  a  variety  of  possible 
wavelengths  ranging  from  4  to  30  m. 


4.  DISCUSSION  OF  SOURCE 
CHARACTERISTICS 

To  understand  the  direction  estimates  of  mov¬ 
ing  seismic-acoustic  sources  it  is  instructive  to 
first  consider  the  spectral  properties  of  simpler 
sources.  We  therefore  begin  by  looking  at  impul¬ 
sive  acoustic  and  seismic  wave  fields  generated  by 
a  .45  caliber  (blank)  pistol  shot  and  a  vertically 
delivered  sledgehammer  blow.  The  pistol  shot 
was  recorded  around  1500  on  16  February,  and  the 
hammer  blow  around  1730  on  the  same  day. 

All  estimates  of  the  phase  velocities  discussed 
in  this  section  are  based  on  a  time-domain  move- 
out  method.  Each  waveform  in  the  time-domain 
signal  vector  is  synchronized  by  specifying  the 
source  direction  (0)  and  phase  velocity  (V0)  of  the 
waveforms.  Since  the  spatial  distribution  of  the 
array  elements  is  known,  we  can  calculate  the 
relative  differences  in  arrival  times  between  all 
elements  in  the  array.  The  time  shift  may  be  ex¬ 
pressed  as 
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Figure  16.  Microphone  subarray  response  to  a  .45  caliber  blank  pistol  shot. 


A tm  =  ^X™+  - 6)  ,  (32) 

where  0Cm  is  the  angular  location  of  the  mth  sensor. 
Given  a  digital  record  of  an  ideal  plane  wave 
source,  the  accuracy  of  eq  32  will  be  limited  by  the 
time-domain  sample  interval. 

4.1.  Impulsive  acoustic  source 

The  first  field  example  discussed  is  the  .45  cali¬ 
ber  pistol  firing  blanks  toward  the  array  from  the 
closest  point  of  approach  (CPA).  The  time-domain 
microphone  and  geophone  response  to  a  single 
pistol  shot  are  shown  in  Figures  16  and  17.  The 
microphone  waveforms  have  phase  velocities  of 
approximately  330  m  s-1.  The  seismic  waveforms 


(Figure  17)  have  two  distinct  arrival  phases.  The 
first  arrival  is  a  roughly  75-Hz  waveform  that  has 
a  phase  velocity  of  approximately  330  m  s-1.  It  has 
the  distinct  appearance  of  the  acoustic  impulse 
observed  in  the  microphone  waveforms.  Later  in 
the  geophone  time  series,  the  waveforms  shift  to 
25  Hz  and  have  a  phase  velocity  of  approximately 
220  ms-1. 

The  frequency  domain  spectra  for  these  signal 
vectors  are  given  in  Figures  18  and  19.  Each  spec¬ 
trum  is  based  on  a  2000-point  series  sampled  at 
0.0005  s,  segmented  into  3  blocks  of  1024  points, 
and  overlapped  by  53%.  A  Blackman  window 
taper  was  applied  to  each  block. 

The  primary  features  to  note  in  the  microphone 
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Figure  17.  Time-domain  vertical  geophone  subarray  response  to  a  .45  caliber 
blank  pistol  shot. 
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Figure  18.  Microphone  spectrum 
for  .45  caliber  pistol  shot:  (left)  Lin¬ 
ear  scaling  and  (right)  logarithmic 
scaling.  The  large  spectral  depres¬ 
sion  at  60  Hz  is  a  result  of  notch 
filtering  during  the  recording  pro¬ 
cess.  The  180-Hzpeak  is  probably  a 
harmonic. 


spectra  are  the  broadband  nature  of  the  signal  and 
the  presence  of  two  high-amplitude  bins  at  45  and 
70  Hz,  separated  by  a  deep  spectral  low  around  60 
Hz  (produced  from  a  notch  filter).  The  appearance 
of  the  geophone  response  is  considerably  more 


complex  and  may  be  a  result  of  several  propaga¬ 
tion  factors  such  as  modal  effects  associated  with 
surface  wave  propagation  and  acoustic-to-seismic 
coupling.  The  majority  of  the  energy  in  vertical- 
geophone  spectra  is  contained  in  a  band  between 


Figure  19.  Spectra  for  three  center 
elements  of  the  geophone  subarray 
for  a  .45  caliber  pistol  shot. 
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Figure  20.  Vertical  geophone  subarray  response  to  a  stack  of  eight 
hammer  blows .  The  element  number  and  peak-to-peak  amplitude  are 
listed  to  the  left  of  each  trace.  The  amplitudes  are  arbitrarily  normalized 
to  the  trace  from  element  9. 


10  and  30  Hz,  with  an  additional  spectral  peak 
near  70  Hz. 

4.2.  Impulsive  seismic  source 

Figure  20  shows  a  series  of  vertical-component 
geophone  waveforms  obtained  by  digitally  stack¬ 
ing  eight  vertically  delivered  sledgehammer  blows. 
The  source  position  is  located  on  the  road  at  a  point 
50  m  to  the  southwest  of  the  CPA  at  a  direction  of 
-78°.  Using  eq  32  it  is  found  that  the  phase  velocity 
is  approximately  220  m  s-1.  The  power  spectra  for 
the  three  center  elements  from  this  array  (Fig.  21) 
indicate  that  the  signal  is  broadband  with  peaks  at 
20  and  40  Hz. 

The  pistol  and  sledgehammer  sources  indicate 
that  the  dominant  seismic  energy  carrier  is  a  sur¬ 
face  wave  propagating  with  a  phase  velocity  of 
approximately  220  m  s-1  and  having  a  mode  fre¬ 
quency  of  between  20  and  30  Hz.  Analysis  of  the 
pistol  shot  records  clearly  shows  that  high-energy 
acoustic  sources  couple  strongly  with  the  ground 
to  produce  surface  waves.  This  is  consistent  with 
the  previous  findings  of  Albert  (1989).  In  addition 
to  the  acoustically  induced  seismic  surface  wave, 
the  geophone  is  also  excited  by  the  airborne  pres¬ 
sure  pulse.  These  propagation  modes  can  interfere 
and  may  be  a  source  of  difficulty  for  beamforming 
processors. 


4.3.  Moving  seismic-acoustic  sources 

Recordings  of  a  tank  were  also  made  as  a  U.S. 
Army  M60  tank  was  driven  up  and  down  the  test 
track  at  various  speeds.  A  malfunction  with  the 
recording  equipment  prevented  continuous  sig¬ 
natures  from  being  obtained,  so  a  series  of  1-s-long 
signature  snapshots  was  recorded.  The  10  tank 
signatures  examined  here  were  recorded  between 
1256  and  1259  on  17  February. 

We  now  discuss  the  spectral  properties  of  an 
M60  tank  moving  at  approximately  4.5  m  s-1  (10 
mph).  The  vehicle's  direction  (0)  during  this  record 
interval  was  approximately  260°.  Figures  22  and 
23  present  the  time-domain  waveforms  observed 
by  the  vertical-geophone  and  microphone 
subarrays,  respectively.  In  the  vertical-geophone 
waveforms  there  is  a  high  degree  of  both  temporal 
and  spatial  variation  in  the  signals.  The  peak-to- 
peak  amplitude  variation  across  the  array  ranges 
up  to  6.6  dB.  Since  the  spatial  variations  between 
array  elements  were  not  consistent  from  record  to 
record,  it  is  unlikely  that  the  variations  are  due  to 
geophone  coupling  to  the  ground  or  to  variation  in 
the  snow  cover. 

Spectra  were  estimated  using  a  total  signal 
length  of  2000  points  segmented  into  three  blocks 
of  1024  points.  Each  block  is  overlapped  by  53% 
and  tapered  using  a  Blackman  window.  Figure  24 
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Figure  21.  Hammer  blow  spectra 
from  center  array  elements  for  ver¬ 
tical  geophone  subarray. 
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Figure  22.  Time-domain  vertical 
geophone  subarray  response  to 
M60  tank  moving  at  4.5  m  sr1. 
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Figure  23.  Time-domain  micro¬ 
phone  subarray  response  to  M60 
tank  moving  at  4.5  m  s~l. 
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Figure  24.  Spectra  of  two 
center  elements  of  vertical 
geophone  subarray  for  an 
M60  tank  moving  at  4.5  m 
S“2.  Spectra  are  strongly 
peaked  at  11  and  29  Hz. 


Frequency  (Hz) 


Frequency  (Hz) 


shows  the  vertical-geophone  power  spectra  of  el¬ 
ements  1  and  2  from  Figure  22.  The  spectra  are 
strongly  peaked  at  11  and  29  Hz.  Plane-wave 
modeling  at  these  two  frequencies  with  a  phase 
velocity  of  220  m  s-1  and  the  Grayling  array  geom¬ 
etry  indicates  that  signal  interference  may  be  the 
cause  of  the  amplitude  variation  across  the  array. 


Interference  effects  may  also  be  present  between 
acoustic  energy  coupled  near  to  the  geophone  and 
surface  waves  excited  near  the  source.  This  hy¬ 
pothesis  is  partially  supported  by  the  presence  of 
a  strong  spectral  peak  at  29  Hz  in  both  the  micro¬ 
phone  and  vertical-geophone  spectra  (Figures  24 
and  25). 


Figure  25.  Spectra  of  micro¬ 
phone  subarray  response  to 
an  M60  tank  moving  at  4.5 
m  s~l.  Spectra  are  strongly 
peaked  around  3  Hz  and  29 
Hz. 
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5.  PROCESSING  RESULTS  FOR 
FIELD  RECORDS 

This  section  presents  ML  wavenumber  esti¬ 
mates  of  Grayling  field  data  for  the  .45  caliber 
blank  pistol,  the  sledgehammer,  and  the  M60  tank 
moving  at  4.5  m  s”1.  All  wavenumber  estimates 
were  formed  over  a  regular  grid  space  containing 
2601  points.  All  spectral  estimates  contain  various 
degrees  of  bias  and  variance.  To  demonstrate  the 
degree  to  which  frequency  domain  bias  and  vari¬ 
ance  can  affect  wavenumber  estimates,  we  show 
several  results  that  use  different  OBAFFT  param¬ 
eters  on  the  same  M60  signal  vector. 

5.1.  Criteria  used  to  select  processed  frequencies 

Several  criteria  were  used  to  select  the  frequency 
bin  at  which  a  wavenumber  estimate  is  formed. 
The  most  important  of  these  processing  criteria 
was  the  average  array  coherence  spectra,  which  is 
formed  by  calculating  and  summing  the  magni¬ 
tude-squared  coherence  spectra  (rave)  between  all 
possible  channel  combinations  and  then  normal¬ 
izing  by  the  number  of  combinations.  This  opera¬ 
tion  is  given  as 

M  |Xj*(/)»Xj*< (fl| 
j=i+l 

^  2  nblksM- 1 

Tave  =  (nblks)(M2-M)  h  h  (33) 

In  general,  wavenumber  spectra  which  were 
processed  at  frequency  bins  with  coherences 
greater  than  0.65  and  a  coherence  variance  across 
immediately  adjacent  bins  of  less  than  0.025,  pro¬ 
duced  reasonable  wavenumber  estimates.  How¬ 
ever,  it  was  found  that  wavenumber  estimates 
could  still  vary  by  a  considerable  degree  even  with 
average  coherences  of  0.9  (usually  this  was  accom¬ 
panied  by  high  coherence  variance  across  adjacent 
bins) .  Such  variation  can  be  explained  by  recogniz¬ 
ing  that  the  coherence  spectrum  is  also  susceptible 
to  bias  and  variance  errors. 

The  average  array  coherence  spectrum  is 
strongly  affected  by  the  OBAFFT  parameters.  The 
Blackman  taper  normally  produced  a  smoothly 
varying  coherence  spectrum  with  peaks  ranging 
over  several  frequency  bins.  In  contrast,  a  boxcar 
taper  often  produced  an  erratic  coherence  spec¬ 
trum.  In  general,  peaks  in  the  coherence  spectra 
greater  than  0.7,  which  varied  smoothly  through 
several  bins,  produced  acceptable  wavenumber 


estimates.  Wavenumber  estimates  based  on  seis¬ 
mic  signals  from  the  moving  tracked  vehicles 
tended  to  be  reliable  at  frequencies  between  9  and 
15  Hz  despite  the  occurrence  of  peak  spectral 
energies  in  the  vicinity  of  30  Hz.  The  30-Hz  bins 
rarely  produced  satisfactory  wavenumber  esti¬ 
mates. 

Several  factors  were  useful  in  qualitatively  de¬ 
termining  the  accuracy  of  wavenumber  spectra. 
The  most  important  was  the  presence  of  a  strong 
array  factor.  This  was  particularly  relevant  when 
assessing  the  quality  of  wavenumber  estimates 
based  on  small  numbers  of  OBAFFT  blocks.  An¬ 
other  key  indicator  proved  to  be  the  estimate  of  the 
phase  velocity.  An  estimated  phase  velocity  con¬ 
sistent  with  the  surface  wave  velocity  found  by  the 
time-domain  move-out  analysis  was  usually  a  good 
indicator  of  an  accurate  wavenumber  estimate.  In 
addition,  if  the  spatial  correlation  is  normalized 
according  to  eq  25  then,  when  using  the  Bartlett 
method,  a  peak  power  value  of  0.6  or  greater  often 
accompanied  an  accurate  estimate. 

5.2.  Impulsive  sources 

The  first  series  of  wavenumber  estimates  is 
based  on  the  vertical-geophone  subarray  response 
to  the  .45  caliber  blank  fired  from  the  CPA  at  0  = 
-40°.  The  waveforms  of  this  signal  vector  are  shown 
in  Figure  17.  Figure  26  shows  the  average  array 
amplitude  and  coherence  spectra,  based  on  a  2000- 
point  series  with  three  blocks.  Each  block  contains 
1024  points,  is  overlapped  by  53%,  and  is  tapered 
using  a  Blackman  window 

Recall  from  section  4  that  there  are  two  distinct 
arrivals  in  this  record.  The  wavenumber  estimate 
shown  in  Figure  27  is  the  Capon  maximum-likeli¬ 
hood  estimate  processed  at  78  Hz.  This  frequency 
was  chosen  since  it  produces  a  wavelength  com¬ 
parable  to  the  spacing  of  the  center  elements  of  the 
geophone  array  with  an  acoustic  propagation  ve¬ 
locity.  The  time-domain  interval  considered  was 
between  0.1  and  0.25  s.  This  corresponds  to  a  time 
window  that  brackets  the  initial  impulsive  arrival 
seen  in  Figure  17.  There  were  230  sample  points  in 
this  interval.  The  estimate  was  formed  using  a 
block  length  of  128  points  with  each  block  over¬ 
lapped  60%,  giving  a  total  of  three  blocks.  Boxcar 
windowing  was  applied.  The  coherence  at  78  Hz 
was  0.7.  The  wavenumber  spectrum  (Fig.  27)  pro¬ 
duces  an  estimated  source  direction  of  -40°  and  a 
phase  velocity  of  338  m  s_1.  This  is  consistent  with 
the  direction  of  the  pistol  shot  and  the  acoustic 
propagation  velocity.  Note  that  the  convention 
followed  plots  the  true  vector  orientation  of  the 
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Figure  26.  Average  array  spectra  for  vertical 
geophone  subarray  to  .45  caliber  pistol  shot. 


Figure  27.  ML  wavenumber  spec¬ 
trum  at  78  Hz  for  vertical  geophone 
subarray.  The  location  of  the  spectral 
maxima  indicates  a  source  direction 
of  0  =  -40°  and  a  phase  velocity  of 
338  ms-l 


East-West  Component  of  Wavenumber  (nr1) 


Linear  grey  scale:  White  Regions  =  -3  dB 
Black  Regions  =  -42  dB 
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Figure  28.  Estimated  ML  wave- 
number  spectrum  at  23  Hz  for  verti¬ 
cal  geophone  subarray.  Location  of 
spectral  peak  gives  a  source  direction 
of  6  -  -37°  and  a  phase  velocity  of 
223  m  s"1. 


Linear  grey  scale:  White  Regions  =  -3  dB 
Black  Regions  =  -42  dB 


wavenumber.  Thus,  all  incoming  source  wave- 
numbers  are  oriented  180°  away  from  the  source 
direction. 

Figure  28  is  a  23-Hz  ML  wavenumber  estimate 
of  the  same  pistol  shot  record.  The  time  interval 
considered  encompasses  the  later  portion  of  the 
wave  train  between  0.2  and  0.5  s.  There  were  500 
points  in  this  interval.  The  time-domain  move-out 
analysis  indicated  a  phase  velocity  of  roughly 
220  m  s'1  in  this  interval.  A  block  length  of  256 
points  overlapped  by  60%  was  used,  giving  four 
blocks  in  the  estimated  spatial  correlation  matrix. 
A  boxcar  window  taper  was  applied.  The  average 
array  coherence  in  the  23-Hz  bin  was  0.8.  The 
wavenumber  spectra  in  Figure  28  give  a  source 
direction  of  0  =  -37°  and  a  phase  velocity  of  223  m 
s_1.  These  values  agree  with  the  known  source 
direction  and  the  time-domain  estimates  of  the 
propagation  velocity. 

The  last  beamformer  response  for  the  pistol 
shot  is  shown  in  Figure  29.  This  is  a  broadband  ML 
estimate  processed  on  an  800-point  time  interval 
between  0.1  and  0.5  s.  A  block  length  of  128  points 
was  chosen  with  no  block  tapering.  The  block 
overlap  was  40%,  giving  nine  blocks.  The  fre¬ 
quency  band  evaluated  was  from  15  to  78  Hz. 


Selection  of  these  parameters  allowed  the  spectral 
estimate  of  each  block  to  retain  a  high  degree  of 
independence  with  emphasis  on  the  resolution  of 
specific  frequency  bins.  In  Figure  29  there  are  three 
spectral  maxima  that  have  energy  levels  greater 
than  -3  dB.  All  three  trend  along  a  high-amplitude 
"ridge/'  indicating  a  source  direction  of  roughly 
-40°.  The  wavenumber  peak  with  the  smallest 
magnitude  yields  a  phase  velocity  of  223  m  s-1  if 
normalized  by  a  frequency  of  23  Hz.  The 
wavenumber  peak  with  the  intermediate  magni¬ 
tude  yields  a  phase  velocity  of  342  m  s”1  when 
normalized  by  78  Hz. 

The  above  series  of  wavenumber  spectra  and 
the  time-domain  plots  clearly  show  the  existence 
of  an  acoustically  induced  seismic  surface  wave. 
In  considering  the  ML  wavenumber  spectra  it 
should  be  remembered  that  the  wavefield  was 
generated  by  a  nearly  ideal  stationary  point  source. 
The  OBAFFT  parameters  used  in  the  ML 
wavenumber  spectra  were  chosen  through  trial 
and  error,  and  in  many  circumstances  reasonable 
parameters  did  not  yield  intelligible  beam  re¬ 
sponses.  The  above  set  of  wavenumber  estimates 
(Fig.  27  to  29)  serve  as  a  testament  to  the  impact  of 
the  frequency  domain  estimation  methodology  on 
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Figure  29.  Broadband  ML  wave- 
number  spectrum  for  vertical  geo¬ 
phone  subarray.  The  processed  fre¬ 
quency  band  is  15  to  78  Hz.  High- 
amplitude  linear  ridge  along  a 
wavenumber  space  direction  ofl  40  ° 
indicates  a  source  position  of -40°. 
The  locations  of  spectral  maxima 
indicate  waves  propagating  at  223 
and  342  m  srl. 


Linear  grey  scale:  White  Regions  =  -3  dB 
Black  Regions  =  -42  dB 


the  resulting  wavenumber  spectra. 

The  next  ML  wavenumber  spectra  we  consider 
are  for  the  vertical-geophone  signal  vector  pro¬ 
duced  by  a  vertically  delivered  sledgehammer 
blow  (see  Fig.  20).  The  source  location  was  on  the 
road  50  m  south  of  the  CPA  at  0  =  -78°.  Figure  30 
displays  the  average  array  power  and  coherence 
spectra.  The  OBAFFT  parameters  used  for  the 
power,  coherence,  and  wavenumber  spectra  were 
a  signal  length  of 2000  points  segmented  into  three 
blocks  of  1024  points,  each  overlapped  by  53%  and 
tapered  with  a  Blackman  window.  The  95%  vari¬ 
ance  bounds  were  R  •  (0.46, 3.96)  with  5.51  d.f.  The 
25-Hz  ML  wavenumber  spectra  based  on  this 
record  are  shown  in  Figure  32.  The  average  array 
coherence  at  25  Hz  is  0.64.  The  location  of  the 
spectral  maximum  in  the  wavenumber  estimate 
gives  0  =  -74°  with  a  phase  velocity  of  232  m  s“L 

5.3.  Moving  sources 

The  first  source  we  consider  (Figure  22)  is  an 
M60  tank  travelling  along  the  road  from  south¬ 
west  to  northeast  at  4.5  m  s-1.  This  was  recorded 
when  the  vehicle  was  300  m  to  the  southwest  of  the 
CPA.  The  OBAFFT  parameters  used  to  produce 
the  wavenumber,  power,  and  coherence  spectra 


Frequency  (Hz) 

Figure  30.  Average  array  spectra  for  vertical  geophone 
subarray  to  stack  of  eight  hammer  blows. 
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Figure  31.  Estimated  ML  wave- 
number  spectrum  at  25  Hz  for  ver¬ 
tical  geophone  subarray  for  ham¬ 
mer  blows.  Peak  amplitude  gives  a 
source  direction  of  6  =  -74° and  a 
phase  velocity  of 232  m  s~L 


Linear  grey  scale:  White  Regions  =  -3  dB 
Black  Regions  =  -42  dB 


are  a  series  length  of  2000  points  segmented  into 
three  blocks  of  1024  points,  each  overlapped  by 
53%,  and  tapered  with  a  Blackman  window  func¬ 
tion.  The  95%  confidence  bounds  are  R  •  (0.46, 3.96) 
with  5.51  d.f.  Figure  32  shows  the  average  array 
power  and  coherence  spectra.  There  is  a  peak  in 
the  power  spectrum  at  11  Hz  that  has  a  coherence 
of  0.94.  The  11-Hz  ML  wavenumber  spectrum 
(Fig.  33)  estimates  the  source  direction  as  0  =  258° 
with  a  phase  velocity  of  219  m  S”1. 

The  result  shown  in  Figure  33  is  typical  of  ML 
wavenumber  estimates  obtained  from  the  M60  in 
the  9-  to  15-Hz  band  using  the  vertical-geophone 
subarray.  To  demonstrate  this,  consider  Figure  34, 
which  is  a  composite  of  direction  and  velocity 
magnitudes  estimated  from  ML  wavenumber  spec¬ 
tra  for  a  consecutive  series  of  M60  vehicle  records. 
The  location  of  each  direction  position  is  sequen¬ 
tially  labeled  1  to  16.  All  records  in  this  composite 
are  from  the  same  M60  drive-by  series,  with  the 
tank  moving  from  southwest  to  northeast  along 
the  road  at  a  constant  speed  of  4.5  m  s-1.  The  most 
distant  records  in  this  series  were  at  positions  1 
and  16,  which  correspond  to  vehicle  ranges  of  500 
and  300  m,  respectively.  At  position  11  the  tank  is 
50  m  from  the  array's  coordinate  origin. 


Figure  32.  Average  array  spectra  for  vertical  geophone 
subarray  to  M60  tank  moving  at  4.5  m  s-L 
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Figure  33.  Example  of  an  ML 
wavenumber  estimate  from  a  verti¬ 
cal  geophone  array  at  11  Hz  for  an 
M60  tank  moving  at  4.5  m  s~J. 
Coherence  was  0.94,  the  spectral 
maximum  yields  a  velocity  of  219  m 
s -1,  and  the  source  direction  =  258  °. 


East-West  Component  of  Wavenumber  (rrr1) 


Linear  grey  scale:  White  Regions  =  -3  dB 
Black  Regions  =  -42  dB 
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Figure  34.  Composite  plot  of  estimated  source  direction  and  velocity 
magnitudes  for  a  series  ofM60  tank  records  observed  by  the  vertical 
geophone  subarray . 
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The  composite  uses  a  range  of  temporal  fre¬ 
quencies  from  9.77  to  29  Hz  with  the  majority  of 
the  records  using  the  11-Hz  bin.  All  the  estimates 
used  2000-point  signal  lengths,  segmented  into 
three  blocks  of  1024  points,  each  overlapping  the 
pluvious  by  53%.  A  Blackman  taper  was  used  for 
10  of  the  wavenumber  estimates,  with  the  remain¬ 
ing  five  estimates  using  a  boxcar  window  taper. 
Using  the  boxcar  window,  the  95%  confidence 
bounds  on  the  crosspower  matrix  were 
R  •(0.42,5.85)  with  3.9  d.f.  The  decision  to  use 
either  a  boxcar  or  Blackman  taper  was  based  on 
the  coherence  spectra  and  the  consistency  of  the 
velocity  estimates  of  the  resulting  wavenumber 
spectra. 

5.4.  Affect  of  spectral  estimation  parameters  on 
beam  response 

In  the  preceding  sections  we  demonstrated  that 
ML  wavenumber  spectra  from  small  seismic  ar¬ 
rays  can  produce  accurate  direction  estimates  with 
careful  choice  of  transformation  parameters.  An¬ 
other  result  is  that  small  differences  in  the  fre¬ 
quency-domain  processing  parameters  often  have 
large  effects  on  the  wavenumber  estimate.  Such 
effects  are  termed  a  wavenumber  bias.  In  this 
section  several  examples  of  the  wavenumber  bias 
phenomenon  are  demonstrated  for  the  M60  tank. 


When  using  reasonable  OBAFFT  parameters, 
the  wavenumber  bias  effect  was  negligible  for 
signal  vectors  observed  from  sledgehammer  blows 
or  blank  pistol  shots.  Thus,  the  wavenumber  bias 
problem  is  most  likely  associated  with  the  charac¬ 
teristics  of  the  M60.  Given  the  difficulty  in  produc¬ 
ing  accurate  wavenumber  estimates  from  the  ver- 
tical-geophone  array  in  the  29-Hz  peak  power  bin 
of  the  M60  source,  and  the  coincident  occurrence 
of  large  energy  levels  at  29  Hz  in  the  microphone 
and  geophone  power  spectra  (Fig.  25  and  26),  it  is 
tentatively  proposed  that  the  bias  effect  may  be 
aggravated  by  interference  between  surface  waves 
and  acoustic  energy  coupling  near  the  geophone 
and  source  movement. 

5.5.  Spatial  frequency  bias  from  moving  sources 

To  show  how  changes  in  processing  param¬ 
eters  can  affect  wavenumber  estimation  for  the 
moving  M60,  we  consider  a  series  of  wavenumber 
estimates  using  the  same  data  set  but  applying 
slightly  different  OBAFFT  parameters.  In  each 
case  it  will  be  shown  how  the  estimation  param¬ 
eters  lead  to  biased  wavenumber  estimates  or  to  a 
reduction  in  the  SNR  as  measured  by  the  back¬ 
ground  energy  level  of  the  wavenumber  spectrum. 

Consider  the  signal  vector  shown  in  Figure  35 
for  an  M60  tank  moving  along  the  road  at  4.5 
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Figure  35.  Time-domain  response  of  vertical  geophone  subarray  to  M60  tank 
moving  at  4.5  m  s~2  along  the  road  approximately  100  m  southeast  of  the  CPA. 
These  waveforms  correspond  to  position  8  in  the  composite  plot  shown  in 
Figure  35. 
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m  s*”1  from  southwest  to  northeast.  An  11-Hz  ML 
wavenumber  spectrum  is  shown  in  Figure  36.  It 
was  formed  using  a  series  length  of  2000  points 
segmented  into  20  blocks;  each  block  had  1024 
points,  was  tapered  by  a  Blackman  window,  and 
overlapped  the  preceding  block  by  95%.  The  95% 
confidence  bounds  on  the  variance  of  the 
crosspower  matrix  was  R«  (0.38, 9.39)  with  2.83 
d.f.  The  average  array  coherence  at  11  Hz  was  0.92. 
The  variance  estimate  leads  one  to  expect  that  no 
significant  improvement  in  the  wavenumber  spec¬ 
tra  will  result  from  such  highly  overlapped  blocks. 
This  is  not  the  case,  as  will  be  shown  later. 

The  wavenumber  spectrum  in  Figure  36  is  no¬ 
table  in  several  respects.  It  has  an  extremely  sharp 
spectral  maximum  with  background  energy  levels 
as  low  as  -45  dB.  The  highest  side  lobe  in  wave- 
number  space  is  30  dB  below  the  maximum.  The 
location  of  the  spectral  maximum  gives  a  source 
direction  of  258°  and  a  propagation  velocity  of  219 
ms-1. 

The  ML  wavenumber  estimate  shown  in  Figure 
38  uses  the  same  signal  vector  as  that  of  Figure  36. 
In  this  spectrum  the  OBAFFT  parameters  are  iden¬ 
tical  except  the  blocks  are  boxcar  windowed.  The 
variance  estimate  is  R  •  (0.42,5.83)  with  3.89  d.f.  at 


the  95%  confidence  level.  In  this  extremely  over¬ 
lapped  region  the  variance  bounds  are  slightly 
better  than  those  found  for  the  Blackman  window 
function  used  in  the  previous  example.  Figure  37 
is  a  notable  wavenumber  estimate  in  that  it  has  an 
exceptionally  peaked  spectral  maximum,  which 
indicates  a  source  direction  of  258°  and  a  phase 
velocity  of  219  m  s~L  There  are  several  subtle 
differences  between  Figures  37  and  38.  The  first  is 
the  occurrence  of  a  relatively  pronounced  spectral 
side  lobe  in  Figure  37  that  is  18  dB  smaller  than  the 
spectral  maximum.  In  addition,  the  background 
energy  level  in  Figure  37  is  approximately  -33  dB 
while  that  of  Figure  36  was  approximately  -45  dB. 

The  next  two  wavenumber  estimates  investi¬ 
gate  the  impact  of  reducing  the  number  of  blocks 
in  the  OBAFFT.  In  Figures  38  and  39  we  process  a 
signal  length  of  2000  points,  blocked  into  seg¬ 
ments  of  1024  points,  with  each  overlapped  by 
53%,  giving  a  total  of  three  blocks.  In  the  case  of 
Figure  38  we  apply  a  Blackman  window  taper  to 
each  block,  whereas  the  wave  number  estimate 
shown  in  Figure  39  used  a  boxcar  window  func¬ 
tion. 

In  Figure  38  the  spectral  maximum  is  sharply 
peaked,  indicating  a  source  direction  of  258°  and  a 
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Figure  36.  ML  wavenumber  esti¬ 
mate  at  11  Hz  for  waveforms  shown 
in  Figure  35:  Blackman  taper ,  20 
blocks  (at  95%  overlap ),  segment 
length  =  1024 ;  produces  a  source 
direction  of 258°  with  a  velocity  of 
219  m  s< 


Linear  grey  scale:  White  Regions  =  -3  dB 
Black  Regions  =  -42  dB 
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Figure  37.  ML  wavenumber  esti¬ 
mate  at  11  Hzfor  waveforms  shown 
in  Figure  35:  boxcar  taper,  20  blocks 
(at  95%  overlap),  segment  length  = 
2 024;  produces  a  source  direction  of 
258°  with  a  velocity  of  219  m  S"1. 


East-West  Component  of  Wavenumber  (nr1) 


Linear  grey  scale:  White  Regions  =  -3  dB 
Black  Regions  =  -42  dB 


phase  velocity  of  219  m  s_1.  The  highest  secondary 
side  lobe  is  at  -21  dB  with  an  overall  background 
energy  level  of -36  dB.  In  contrast,  the  wavenumber 
estimate  shown  in  Figure  39  indicates  a  source 
direction  of  255°  and  a  phase  velocity  of  271  m  s_1 . 
The  3°  direction  difference  and  52  m  s-1  velocity 
difference  between  Figures  38  and  39  is  a  mild 
example  of  the  wavenumber  bias  effect.  In  addi¬ 
tion,  applying  the  boxcar  window  effectively  re¬ 
duces  the  SNR  of  the  beamformer  in  the  sense  that 
the  background  energy  level  of  Figure  38  increases 
by  6  dB  in  Figure  39. 

Comparing  Figures  36  and  37  with  38  and  39  we 
note  that  there  seems  to  be  significant  benefit  in 
using  a  larger  number  of  highly  overlapped  blocks. 
Each  estimate  of  the  phase  difference  between  two 
channels  for  a  given  block  of  data  is  calculated 
from  a  comparison  of  two  vectors.  Each  vector  has 
a  magnitude  and  a  phase,  both  of  which  vary  with 
time,  and  each  contains  a  small  amount  of  random 
noise.  The  amplitude  at  a  given  frequency  will 
vary  slowly  across  blocks  for  a  given  channel, 
indicating  highly  correlated  blocks.  The  phase 
angle  variation  across  blocks  changes  (rotational 
orientation)  rapidly  compared  with  the  amplitude 
variation,  indicating  a  lower  degree  of  correlation 


between  blocks.  Thus,  if  one  is  interested  in  esti¬ 
mating  the  phase  angle  difference  (i.e.,  calculat¬ 
ing  a  spatial  correlation  matrix)  between  chan¬ 
nels  in  an  array  and  the  signal  duration  is 
fixed,  it  is  advantageous  to  use  a  large  number  of 
highly  overlapped  blocks  (compare  Fig.  36  and 
Fig.  38). 

Continuing  with  the  same  data,  we  contrast  the 
previous  four  satisfactory  wavenumber  spectra 
with  another  pair  of  11-Hz  ML  wavenumber  esti¬ 
mates  formed  using  a  single  block  of  512  points. 
This  is  theoretically  more  than  twice  the  signal 
duration  needed  to  adequately  observe  an  11-Hz 
waveform  propagating  at  220  m  s-1. 

The  11-Hz  ML  wavenumber  spectrum  in  Fig¬ 
ure  40  applied  a  Blackman  window  taper  to  the 
block  of  data,  and  the  11-Hz  ML  wavenumber 
spectrum  in  Figure  41  used  a  boxcar  window 
function.  Figure  40  has  an  acceptably  peaked  spec¬ 
tral  maximum  indicating  a  source  direction  of 258° 
and  a  phase  velocity  of  219  m  s~L  However,  this 
estimate  has  a  considerably  higher  background 
energy  level  (-18  dB)  compared  with  the  previous 
examples.  The  highest  spectral  side  lobe  has  an 
energy  level  of  -12  dB.  In  a  rather  dramatic  con¬ 
trast,  we  note  that  the  ML  wavenumber  spectrum 


30 


Figure  38.  ML  wavenumber  esti¬ 
mate  at  11  Hz  for  waveforms  shown 
in  Figure  35:  Blackman  taper,  three 
blocks  (at  53%  overlap),  segment 
length  =  1024;  produces  a  source 
direction  of 258°  with  a  velocity  of 
219  m  s^1. 


Figure  39.  ML  wavenumber  esti¬ 
mate  at  11  Hz  for  waveforms  shown 
in  Figure  35:  boxcar  taper,  three 
blocks  (at  53%  overlap),  segment 
length  -  1024;  produces  a  source 
direction  of 255°  with  a  velocity  of 
271  m  srl. 
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Figure  40.  ML  wavenumber  esti¬ 
mate  at  11  Hz  for  waveforms 
shown  in  Figure  35:  Blackman 
taper ,  one  block ,  segment  length  = 
512.  Spectra  give  a  source  direc¬ 
tion  of  25°  and  velocity  of  219 
m  s-1. 


Figure  41.  ML  wavenumber  esti¬ 
mate  at  11  Hz  for  waveforms 
shown  in  Figure  35:  boxcar  taper , 
one  block,  segment  length  =  512. 
Beam  is  unfocused. 


Linear  grey  scale:  White  Regions  =  -3  dB 
Black  Regions  =  -42  dB 


in  Figure  41  is  essentially  unfocused  for  the  boxcar 
taper. 

Figure  41  is  an  extreme  example  of  beamformer 
failure;  however,  it  was  not  uncommon  to  encoun¬ 
ter  similar  problems  in  less  coherent  frequency 
bins  with  more  reasonable  OBAFFT  parameters. 
The  failure  of  the  beamformer  to  focus  is  most 
likely  due  to  the  highly  biased  nature  of  the  fre¬ 
quency-domain  phase  angle  caused  by  energy 
leakage  from  neighboring  frequency  bins.  Such 
bias  effects  are  particularly  problematic  with  the 
boxcar  window.  The  above  examples  show  the 
extent  to  which  the  introduction  of  a  high  degree 
of  bias  and  variance  in  the  frequency  domain 
estimate  impacts  the  wave  number  estimate. 

Considering  the  wavenumber  results  in  this 
section  and  those  given  in  other  sections,  we  make 
several  generalizations  concerning  the  impact  of 
the  OBAFFT  parameters  on  ML  beamformer.  In 
virtually  all  the  examples  discussed,  the  applica¬ 
tion  of  a  strongly  tapered,  low-resolution,  low- 
bias,  time-domain  window  function  (such  as  the 
Blackman  taper)  produced  a  higher  resolution 
wavenumber  estimate  with  lower  background 
energy  levels  when  compared  with  wavenumber 
estimates  using  high-resolution,  high-bias,  time- 
domain  tapers  (such  as  the  boxcar  function). 

Specifying  large  numbers  of  blocks  in  the 
OBAFFT  estimation  of  the  spatial  correlation  ma¬ 
trix  mitigated  the  impact  of  the  high-bias  window 
functions.  Improvements  in  beamformer  perfor¬ 
mance  were  observed  even  when  blocks  were 
overlapped  by  as  much  as  95%.  The  unexpected 
benefit  of  using  highly  overlapped  blocks  in  the 
estimate  of  the  spatial  correlation  matrix  may  be 
explained  by  noting  that  the  phase  angle  for  a 
given  frequency  bin  in  each  block  varies  rapidly 
between  blocks  and  therefore  has  a  much  smaller 
degree  of  correlation  between  blocks  as  compared 
with  power  spectra  estimates.  The  variance  analy¬ 
sis  was  helpful  in  partially  explaining  the 
beamformer's  dependence  on  the  OBAFFT  pa¬ 
rameters.  A  more  complete  explanation  will  have 
to  consider  the  phase  angle  bias  effects  as  well. 

6.  CONCLUSIONS 

Capon's  minimum  variance  (maximum  likeli¬ 
hood)  beamformer  produces  reliable  direction  es¬ 
timates  for  a  U.S.  Army  M60  tank  moving  at  4.5 
m  s_1  using  an  array  of  vertical-component  geo¬ 
phone  signals.  Accurate  direction  estimates  were 
obtained  at  ranges  from  50  to  500  m  from  the  array 


origin.  The  vertical-geophone  responses  to  the 
M60  tank  were  dominated  by  the  presence  of 
Rayleigh  surface  wave  energy  with  spectral  peaks 
at  11  and  30  Hz  and  a  propagation  velocity  of 
220  m  s-1.  Most  wavenumber  estimates  used  the 
highly  coherent  secondary  maxima  in  the  vicinity 
of  the  11-Hz  frequency  bins.  It  was  difficult  to 
obtain  reliable  wavenumber  estimates  in  the  vicin¬ 
ity  of  the  30-Hz  spectral  maxima.  It  is  hypoth¬ 
esized  that  the  low  coherence  in  the  vicinity  of  30 
Hz  is  due  to  interference  at  the  sensor  between 
acoustic  and  seismic  waves. 

Maximum-likelihood  wavenumber  estimates 
of  .45  caliber  blank  pistol  shots  exhibited  two 
distinct  propagation  modes.  The  first  wave  phase 
arrives  with  a  propagation  velocity  of  338  m  s”1 
and  is  interpreted  as  an  acoustic,  excitation  of  the 
ground  in  the  vicinity  of  the  geophone.  The  second 
wave  phase  arrives  with  a  propagation  velocity  of 
approximately  220  m  s-1  and  is  interpreted  as  a 
Rayleigh  surface  wave  that  is  excited  by  acoustic- 
to-seismic  coupling  near  the  source.  Wavenumber 
estimates  of  vertically  delivered  sledgehammer 
blows  also  exhibited  surface  wave  propagation 
modes  that  had  the  same  phase  velocity. 

In  general,  direction  and  velocity  estimates  ob¬ 
tained  from  wavenumber  spectra  showed  good 
agreement  with  known  source  positions  and  ve¬ 
locity  estimates  produced  by  time-domain  move- 
out  analyses.  However,  it  was  observed  that 
changes  in  the  parameters  used  to  estimate  the 
spatial  correlation  matrix  could  produce  wave- 
number  bias  effects.  These  spatial  frequency  bias 
effects  were  strongly  dependent  on  the  type  of 
window  taper  applied  to  the  time  series  data  and 
the  number  of  blocks  used  in  the  estimate  of  the 
spatial  correlation  matrix.  It  was  observed  that 
low-resolution,  low-bias  window  types  produced 
the  highest  resolution  wavenumber  spectra  with 
the  least  wavenumber  bias  and  highest  signal-to- 
noise  ratios.  Utilizing  large  numbers  of  blocks  in 
the  estimate  of  spatial  correlation  matrix  miti¬ 
gated  the  impact  of  window  taper  choice.  In  some 
cases  improvements  were  found  in  the  wave-num¬ 
ber  spectra  even  when  block  overlaps  were  as  high 
as  95%. 
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